IS_CONTROLLER = True

if IS_CONTROLLER: from chimera.controllers import database
else: import database
import chimera.core.log
from chimera.core.event import event
from chimera.util.position import Position
from chimera.util.sextractor import SExtractor
from chimera.util.coord import Coord
from datetime import datetime
from optparse import OptionParser
import ephem
import math
import os
import pyfits
from types import StringType
from chimera.core.chimeraobject import ChimeraObject
from chimera.core.remoteobject import RemoteObject
import copy
import sys
import re
from pyraf import iraf

# These two variables are used globally later on. This is bad style, but it works. FIXME
latitude_from_config = longitude_from_config = 0

class Star:
    def __init__(self, ra = None, dec = None, objname = None):
        self._objname = objname
        self._ra = ra
        self._dec = dec

    def getObjName(self):
        return self._objname

    def setObjName(self, objname):
        self._objname = objname

    def getRA(self):
        return self._ra

    def getDec(self):
        return self._dec

    def display(self):
        if self.getObjName() != None:
            print self.getObjName()
        print "RA:  ", self.getRA()
        print "Dec: ", self.getDec()


class Exposure(RemoteObject):
    #The Exposure class contains information about an entire image.
    def __init__(self, filename=None):
	RemoteObject.__init__(self)
        self._expstars = []
	# Read all the needed information from the FITS header.  	
	if filename != None:
            self._filename = filename
            hdulist = pyfits.open(filename)
	    print '\n'+80*'~'
	    print 'Running photometry on the following file:'
	    self.display()
	    print '\n'
	# Because Epoch is deprecated, we change it to Equinox if used. We also change Epoch OB to Equinox. If there is no Equinox, Epoch or Epoch OB specified we default to the present month.
	currentdate = datetime.now()
	currentequinox = currentdate.year+(currentdate.month/12)
	try:
		hdulist[0].header.update('Equinox',hdulist[0].header['Equinox'])
	except:
		try:
		   hdulist[0].header.update('Equinox',hdulist[0].header['Epoch'])
		   if not IS_CONTROLLER:
		   	print 'Substituted Equinox for Epoch; Epoch is depracated.\n'
		except:
		   try:
		      hdulist[0].header.update('Equinox',hdulist[0].header['Epoch OB'])
		      if not IS_CONTROLLER:
		   	print "Substituted Equinox for Epoch OB; Epoch OB is invalid.\n"
		   except:
		      hdulist[0].header.update('Equinox',currentequinox)
		      if not IS_CONTROLLER:
		   	print "Defaulting to Equinox of "+currentequinox

	# Verify that all keywords are allowed. If there are disallowed keywords, attempt to fix them and warn the user. 
	try:
	    hdulist.verify('silentfix')
	except:
	   print '\nThe following keywords are illegal. Photometry will attempt to continue.\n'
	   hdulist.verify('warn')
	   print '\n'
	
	
        # Get latitude from FITS header or site information
        if 'LATITUDE' in hdulist[0].header:  self._latitude = Coord.fromDMS(hdulist[0].header['LATITUDE'])
        elif 'SITELAT' in hdulist[0].header: self._latitude = Coord.fromD(hdulist[0].header['SITELAT'])
        else:                  
		if not IS_CONTROLLER:
		   	print "Database entry defaulting to supplied image source latitude."
		self._latitude = Coord.fromDMS(latitude_from_config)


        # Get longitude from FITS header or site information
        if 'LONGITUD' in hdulist[0].header:   self._longitude = Coord.fromDMS(hdulist[0].header['LONGITUD'])
        elif 'SITELONG' in hdulist[0].header: self._longitude = Coord.fromD(hdulist[0].header['SITELONG'])
        else:                     
		if not IS_CONTROLLER:
		   	print "Database entry defaulting to supplied image source longitude."		
		self._longitude = Coord.fromDMS(longitude_from_config)

        # Get observation date and time from FITS header	

	try:
		try:
			datestring=str(hdulist[0].header['DATE-OBS'])+'T'+str(hdulist[0].header['UT'])
			try:
				obsdate = datetime.strptime(datestring, '%Y-%m-%dT%H:%M:%S.%f')
			except:
				obsdate = datetime.strptime(datestring, '%Y-%m-%dT%H:%M:%S')
		except:
			note =str(hdulist[0].header['NOTE'])
       			obsdate=self.parse(note)
			try:
		  		obsdate = datetime.strptime(obsdate, '%Y-%m-%dT%H:%M:%S.%f')
			except:
		   		obsdate = datetime.strptime(obsdate, '%Y-%m-%dT%H:%M:%S')
		if not IS_CONTROLLER:
		   		print "Read in observation date successfully."
	except:
		try:        
	   	   try:
			obsdate = datetime.strptime(hdulist[0].header['DATE-OBS'], '%Y-%m-%dT%H:%M:%S.%f')
	   	   except:
			obsdate = datetime.strptime(hdulist[0].header['DATE-OBS'], '%Y-%m-%dT%H:%M:%S')
			if not IS_CONTROLLER:
		   		print "Read in observation date successfully."
        	except:
            	   try:    
			try:
			   obsdate = datetime.strptime(hdulist[0].header['DATE'], '%Y-%m-%dT%H:%M:%S.%f')
			except:
			   obsdate = datetime.strptime(hdulist[0].header['DATE'], '%Y-%m-%dT%H:%M:%S')
	    		
		   	print "\nUsing FITS generation date as observation time. Results may be inaccurate."		
		   	print "Examine fits headers and ensure the presence of a DATE-OBS keyword.\n"
            	   except: 
			obsdate = datetime.now()
		   	
		   	print "\nObservation date or time not found! No correct results possible!"
		   	print "Observation date read unsuccessfully.\n"

	try:
		obs = ephem.Observer()
	except:
		if not IS_CONTROLLER:
			print "\nUnable to create an observer ephemeresis object.\n"
	
	obs.lat = self._latitude.strfcoord('%(d)d:%(m)d:%(s).2f')
	obs.long = self._longitude.strfcoord('%(d)d:%(m)d:%(s).2f')

	# Sets the air pressure to use in refraction calculations near the horizon. Set to local air pressure, or 0 to disable this feature. If deleted calculations will procede using 1010.0 millibars.
	obs.pressure = 0
	if ((obs.pressure==0) and (not IS_CONTROLLER)):
		print '\nNot using PyEphem refraction calculator.\n'
       	obs.date = obsdate

	    #Calculate ST using the sidereal_time method in the python package ephem.

	self._lst = Coord.fromHMS(str(obs.sidereal_time()))
		

	    #All Chimera images have the keyword FILTER in the header.  If you are using images from another observatory, the keyword might be FILTERS	
	try:
		self._filter = hdulist[0].header['FILTER']
	except:
		try:
			self._filter = hdulist[0].header['FILTERS']
		except:
			self._filter = None
	
	# Finally, get the size of the image. This is useful for importing all the stars found in Simbad in the frame.

	self._size = hdulist[0].header['NAXIS1']
	if hdulist[0].header['NAXIS2'] > self._size:
		self._size = hdulist[0].header['NAXIS2']
	
	# Uncomment to print info about the exposure recorded by photometry when photometry is run. Useful in debugging.

	#print '\n\nPhotometry Info:\nFilter: '		
	#print self._filter
	#print'\nLST: '
	#print self._lst
	#print'\nLatitude: '
	#print obs.lat
	#print'\nLongitude: '
	#print obs.long
	#print'\n'

		
    def parse(self, str):
    	if not re.search('\d\d?/\d\d?/\d\d\d\d\ at \d\d?:\d\d?:\d\d?', str): return "Invalid format"
        day = str.partition('/')[2].partition('/')[0]
        month = str.partition(':')[2].partition('/')[0]
        hour = str.partition(' at ')[2].split(':')[0]
        minute = str.partition(' at ')[2].split(':')[1]
        second = str.partition(' at ')[2].split(':')[2]
        if len(month) is 1: month = "0" + month
        if len(day) is 1: day = "0" + day
        if len(hour) is 1: hour = "0" + hour
        if len(minute) is 1: minute = "0" + minute
        if len(second) is 1: second = "0" + second
        return str.partition('/')[2].partition('/')[2].partition(' ')[0] + "-" + month + "-" + day + "T" + hour + ":" + minute + ":" + second

    def getFileName(self):
        return self._filename
    
    def getRA(self):
        return self._ra

    def getDec(self):
        return self._dec

    def getAlt(self):
        return self._alt

    def getAz(self):
        return self._az

    def getLatitude(self):
        return self._latitude

    def getLST(self):
        return self._lst

    def getFilter(self):
        return self._filter

    def appendStar(self, expstar):
        if isinstance(expstar, ExposureStar):
            self._expstars.append(expstar)

    def getAllStars(self):
        return copy.deepcopy(self._expstars)

    def getBrightestStars(self, n = 10):
        added = False
        brightestStars = []
        for expstar in self._expstars:
            i = 0
            added = False
            for stored in brightestStars:
                if expstar.getMagnitude() < stored.getMagnitude():
                    brightestStars.insert(i, expstar)
                    if len(brightestStars) > n:
                        brightestStars.pop()
                        added = True
                    break
                i = i + 1
            if len(brightestStars) < n and not added:
                brightestStars.append(expstar)
        return copy.deepcopy(brightestStars)

    def trimToBrightestStars(self, n = 10):
        self._expstars = self.getBrightestStars(n)
   
    def getSize(self):
    	return self._size

    def display(self):        
	print "File:", self.getFileName()

	
class ExposureStar:
    # The exposure star class contains information about individual stars in an image.  
    def __init__(self, exposure = None, star = None):
	#convert RA/Dec given by SExtractor into Alt/Az.          
	if exposure != None and star != None:

            self._star = star
	    self._exposure=exposure
            raDecTuple = star.getRA(), star.getDec()
            raDec = Position(raDecTuple)
            altAz = Position.raDecToAltAz(raDec, exposure.getLatitude(), exposure.getLST())
            self._alt = altAz.alt
            self._az = altAz.az
            self._Ximage = None
            self._Yimage = None
  

	# Uncomment the following to get more info about each processed star printed.

	#print '\n\n~~~Star Info~~~\nStar: '		
	#print self._star
	#print '\nMore Accurate RA: '
	#print star.getRA()
	#print'\nRA, Dec: '
	#print raDec
	#print'\nAlt, Az: '
	#print altAz
	
    def getStar(self):
    	return self._star


    def getObjName(self):
        return self._star.getObjName()

    def setObjName(self, objname):
        self._star.setObjName(objname)

    def getRA(self):
        return self._star.getRA()

    def getDec(self):
        return self._star.getDec()

    def getAlt(self):
        return self._alt

    def getAz(self):
        return self._az

    def getFlux(self):
        return self._flux

    def setFlux(self, flux):
        self._flux = flux

    def getFWHM(self):
        return self._fwhm
        
    def getXimage(self):
    	return self._Ximage
    
    def getYimage(self):
    	return self._Yimage

    def setFWHM(self, fwhm):
        self._fwhm = fwhm

    def getMagnitude(self):
        return self._mag

    def setMagnitude(self, mag):
        self._mag = mag
        
    def setXimage(self, Ximage):
    	self._Ximage= Ximage
    
    def setYimage(self, Yimage):
    	self._Yimage= Yimage
    	
    def getAirmass(self):
	# Calculation of Airmass based on SExtractor outputs        
	# Zenith angle = PI/2 - Altitude
        # Airmass = sec(Zenith angle)
	
        return (1 / math.cos(math.pi / 2 - self._alt.R))


    def display(self):      
	# Print Photometry results to the terminal.  This is used only if Photometry is being run as a util.	
	if self.getObjName() != None:
            print self.getObjName()
        print "RA:  ", self.getRA()
        print "Dec: ", self.getDec()
        print "Alt: ", self.getAlt()
	print "Az: ", self.getAz()
        print "Mag: ", self.getMagnitude()
        print "Flux:", self.getFlux()
        print "FWHM:", self.getFWHM()
        print "AM:  ", self.getAirmass()
       
	
class Photometry(ChimeraObject):
    exposure = None
    expstars = []
    


    def _getStar(self, ra, dec, objname = None):
	if objname == None:
		return Star(ra, dec)
	else:
        	return Star(ra, dec, objname)

    def run(self, exposure, stars = None, n=30):
	'''Use SExtractor to find the RA/DEC, flux, FWHM, Alt/Az, and magnitude of stars and write the results to the database.  If the list stars is empty, the brightest n (default 30) stars will be written to the database.  To write specific stars to the database, create a file including the ra, dec, and object name of each star.'''
	site_info = self.getManager().getProxy("/Site/0")
	
	global latitude_from_config 	
	latitude_from_config = str(site_info['latitude'])
    	global longitude_from_config
	longitude_from_config = str(site_info['longitude'])
	
        if exposure != None:
            if isinstance(exposure, StringType):
                exposure = Exposure(exposure)
            if isinstance(exposure, Exposure):
                self.exposure = exposure
		
                # Extract point sources and gather data
                sex = SExtractor()
                sex.config['MEMORY_PIXSTACK'] = 9900000 # default 300000
	    	sex.config['BACK_TYPE'] = "AUTO"
	    	sex.config['DETECT_THRESH'] = 3.0
            	sex.config['DETECT_MINAREA'] = 18.0         
	        #sex.config['PHOT_AUTOPARAMS'] = 10,10
                sex.config['PARAMETERS_LIST'] = ['NUMBER',
                                                 'X_IMAGE',     # Pixels
                                                 'Y_IMAGE',     # Pixels
                                                 'X_WORLD',     # D (convert to Coord)
                                                 'Y_WORLD',     # D (convert to Coord)
                                                 'MAG_BEST',    # Magnitude
                                                 'FWHM_IMAGE',  # Pixels (or should it be'FWHM_WORLD'?)
                                                 'FLUX_BEST']   # Count
                sex.run(exposure.getFileName())
                catalog = sex.catalog()
                
                # Load proper IRAF packages and do IRAF photometry.
                
                print "Opening IRAF Packages; ignore the following output:\n"
		iraf.noao()
		iraf.noao.digiphot()
		iraf.noao.apphot()
		xcoords=[]
		ycoords=[]
		magnitudes=[]
		print "\nIRAF Packages successfully opened!"	
		
		
		for item in catalog:
			xcoords.append(item['X_IMAGE'])
			ycoords.append(item['Y_IMAGE'])
		
		# This section contains the code which finds the magnitude of each stars using IRAF, rather than merely Sextractor. It is a major work in progress, especially with Pyraf syntax and data flow. Thus for now it is entirely commented out.
#		coords = self.writeCoordFile(xcoords, ycoords)
#		#self._getIRAFmag(exposure.getFileName())
#		iraf.noao.digiphot.apphot.phot(image=exposure.getFileName(), coords = "coords.coo", centerpars = "", datapars="", interactive = "no")
		#magnitudes = read outputfile
#		self.removeCoordFile()
		#self.removeOutputFile(outputfile)
		
                # Store results
                n = 0
                for item in catalog:
                    # Convert pixel coordinates of every point source to RA/DEC
		    star = self._getStar(Coord.fromD(item['X_WORLD']).toHMS(),
                                         Coord.fromD(item['Y_WORLD']).toDMS()
                                         )
                    try:
                    	a = 1/0#Junk failure try block.
#			magnit = magnitudes[n]
                    except:
                    	magnit = item['MAG_BEST']
                                      
                    n += 1
                    
		    #get photometry info on all stars, if a list of stars was not given	
		    if (stars == None):

			expstar = ExposureStar(exposure, star)
                    	expstar.setFWHM(item['FWHM_IMAGE'])
                    	expstar.setFlux(item['FLUX_BEST'])
                    	expstar.setMagnitude(magnit)
                    	expstar.setXimage(item['X_IMAGE'])
                    	expstar.setYimage(item['Y_IMAGE'])
                    	expstar.setObjName('Star ' + str(item['NUMBER']))
		    	exposure.appendStar(expstar)	
		    else:
			for s in stars:
				# Get photometry on stars with an ra and dec within 2 arcseconds of objects in the list stars		                
				TOLERANCE1 = 0.008
				TOLERANCE2 = 0.0006
				if(abs(Coord.fromHMS(s[0]).D-star.getRA().D)<TOLERANCE1) and (abs(Coord.fromDMS(s[1]).D-star.getDec().D)<TOLERANCE2):
					expstar = ExposureStar(exposure, star)
                    			expstar.setFWHM(item['FWHM_IMAGE'])
                    			expstar.setFlux(item['FLUX_BEST'])
                    			expstar.setMagnitude(magnit)
                    			expstar.setObjName(s[2])
		    			exposure.appendStar(expstar)

		#Trigger photometryCompleted event
		sex.clean(config=True, catalog=False, check=True)

		#Only select the brightest stars if no stars are specified in a list. Default number of stars is 30, but n can be selected when running the script chimera-photo. 
		if stars == None:
			exposure.trimToBrightestStars(n)                	

		daemon = self.getDaemon()
		if daemon == None:		
			self.photometryCompleted(exposure)
			return exposure
		else:
			daemon.connect(exposure)
			self.photometryCompleted(exposure.getProxy())
			return exposure.getProxy()

    # This function represents the core of the processing required to 		
    def getIRAFmag(self, file, skyfile, xcoord, ycoord):
    	return self._getIRAFmag(file, xcoord, ycoord)

    def _getIRAFmag(self, file, xcoord, ycoord):
    	
    	###############################################
    	#
    	# There should probably be some logic in here.
    	#
    	###############################################
    	return None
    	
    # This file takes the x and y coordinate lists from the sextractor catalog file and puts them in an IRAF-readable form.
    def writeCoordFile(self, xlist, ylist):
    	
    	try: os.remove('coords.coo'); print "Removing old coords.coo"
    	except: None
    	try:
    		coords = file('coords.coo','w')
  
    		coords = open('coords.coo', 'w')

    	except:
    		print "Unable to create coordinate file for IRAF photometry!\n"
    	try:
    		n =0
    		while n < len(xlist):
    			coords.writelines(["%s, %s\n" % (xlist[n], ylist[n])])


    			n += 1
    		coords.close()
    	except:
    		print "Unable to write stellar coordinates to coords.coo!\n"
    		
    # This function removes the coordinate file when it is no longer necessary for the sake of not crowding up the system with too many files.	
    def removeCoordFile(self):
    
    	try: os.remove('coords.coo'); print "Getting rid of old coords.coo"
    	except: print "Unable to delete coords.coo, but this should not affect photometry operations."
	
    # This function reads in the output of IRAF's photometry	
    def parseOutputFile(self, outputfile):
    
    	magnitudes = []
    	print "This shouldn't do anything yet."
    	#########################################
    	#
    	# This is where you should include code which takes in the IRAF output file and parses the magnitudes
    	#
    	########################################
    	return magnitudes
    
    # This gets rid of the output file when it is no longer needed. Can be disabled at the cost of having spare files floating around.
    def removeOutputFile(self, outputfile):
    	
    	try: os.remove('SPECIFY OUTPUT FILE')
    	except: print ("\nUnable to remove output file %s, but this should not adversely affect photometry operations." % outputfile)

    @event
    def photometryCompleted(self, exposure):
        pass


